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In this work we have used extensive Monte Carlo simulations and finite size scaling theory to 
study the phase transition in the dipolar Planar Rotator model (dPRM) , also known as dipolar 
XY model. The true long-range character of the dipolar interactions were taken into account by 
using the Ewald summation technique. Our results for the critical exponents does not fit those from 
known universality classes. We observed that the specific heat is apparently non-divergent and the 
critical exponents are v = 1.277(2), /3 = 0.2065(4) and 7 = 2.218(5). The critical temperature was 
found to be Tc = 1.201(1). Our results are clearly distinct from those of a recent Renormalization 
Group study from Maier and Schwabl [PRB 70, 134430 (2004)] and agrees with the results from 
a previous study of the anisotropic Heisenberg model with dipolar interactions in a bilayer system 
using a cut-off in the dipolar interactions [PRB 79, 054404 (2009)]. 
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I. INTRODUCTION 



The Planar Rotator model (PRM) in two dimensions, 
also known as XY model, is known to have a criti- 
cal line in the low temperature region ^, ^. The 
PRM is described by the following Hamiltonian: H = 
-Jj:<.,j>Sr ■ S, = -JE<,,,>(5f^f + SfS]), where 
Si is a two dimensional vector {Sf , Sf) defined in the 
sites i of a two-dimensional lattice and < i,j > means 
that the summation is to be evaluated for nearest neigh- 
bors sites. As a prototype the PRM is expected to 
describe the magnetic properties of ferromagnetic thin 
films where the spins lie in the film plane. Although 
very simple, this model presents some unusual charac- 
teristics, as the absence of spontaneous magnetization 
for any T > 0, which is a consequence of the Mermin 
and Wagner theorem Q . Thus, the system can not have 
a phase transition of order-disorder type, nevertheless, 
there is still a phase transition in the model character- 
ized by a change in the spin-spin correlation function 
behavior. It is observed an algebraic decay of the correla- 
tion function below a characteristic temperature, Tbkt, 
above which the decaying is exponential. Besides that, 
the correlation length is expected to diverge exponen- 
tially as long as Tbkt is approached from above, i.e., 
^(r) - aexp(5/Vr- TBKT) for T > Tbkt, while it 
remains infinity for any T < Tbkt- This transition 
is named Berezinskii-Kosterlitz-Thouless (BKT) phase 
transition ^. Several works, analytical as well as nu- 
merical, dealing with the subject were published since 
the seminal work of Berezinskii and Kosterlitz and Thou- 



* Permanent address: Departamento de Ffsica, Universidade Fed- 
eral de Vir;osa,. .^6570-000 Vigosa, Minas Gerais, Brazil.; Electronic 
address: lucasmol@ufv.bi 



less 1^, H 0, H ^- Besides that, it is also observed 
that the specific heat does not diverge, instead, it has 
a broad maximum at a temperature slightly higher than 
Tbkt 0, [l^, HI, ■ There are two interpretations for 
the mechanism leading to this transition: Berezinskii 
and Kosterlitz and Thouless |^ assume that it is driven 
by a vortex-anti-vortex unbinding mechanism, while Pa- 
trascioiu and Seller |14 were able to obtain the critical 
temperature and predicted the existence of a phase tran- 
sition in the Coulomb gas in any dimension (d > 1) by 
considering that the mechanism responsible for the tran- 
sition is a polymerization of domain walls. (As a matter 
of unification of language we use in this paper the termi- 
nology BKT for this kind of transition). 



However, in order to achieve a deeper insight on the 
magnetic properties of thin films, one has to include dipo- 
lar interactions between the magnetic moments of the lat- 
tice. This inclusion changes the scenario drastically, as 
discussed by Maleev The long-range dipolar inter- 
actions stabilizes the magnetization at low temperatures 
in such a way that an order-disorder phase transition is 
now expected to take place. In a recent paper, Maier and 
Schwabl have analyzed the phase transition in the 
dipolar Planar Rotator model (dPRM) by using renor- 
malization group techniques. Their results indicate that 
the dPRM belongs to a new universality class charac- 
terized by an exponential behavior of the magnetization, 
susceptibility and correlation length. Besides that, the 
specific heat was found to be non-divergent, like occurs 
in the BKT phase transition. In this work, we have used 
extensive Monte Carlo simulations to study the phase 
transition in the dPRM. Our results clearly indicate that 
the transition is of order-disorder type and is character- 
ized by a non-divergent specific heat and unusual critical 
exponents. 
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II. DIPOLAR PLANAR ROTATOR MODEL 
AND MONTE CARLO METHOD 



The model we are interested in consists of a square 
lattice with dimension L x L. At each site we place a 
classical spin variable Si = {Sf,Sf) with 5*^^ — 1. The 
interactions are defined by the following Hamiltonian: 



H = ~J 



<i,]> 



Si ■ Sj 3{Si ■ rij)(Sj ■ Tij) 



(1) 

Here, J > defines a ferromagnetic exchange constant 
and D is the dipolar constant, fij connects sites i and 
J while < i,j > means that the first summation is 
evaluated for nearest neighbors only. For the dipolar 
interactions the summation is evaluated over all pairs 
i ^ j. Periodic boundary conditions have been used in 
the film plane {x and y directions) while open bound- 
ary conditions were ap plie d in the z direction. Ewald 
summation techniques (l7[ have been used to take 
into account the true long-range character of the dipo- 
lar interactions^^. In all simulations we have assumed 
J = 1 and D ~ 0.1 and for these values only ferromag- 
netic configurations were found in the low temperature 
regime. In this work the energy is measured in units of 
JS^ and temperature in units of JS^ /ks- 

Our Monte Carlo procedure consists of a simple 
Metropolis algorithm ^ where one Monte Carlo 
step (MCS) consists of an attempt to assign a new ran- 
dom direction to each spin in the lattice. To equilibrate 
the system we have used 100 x MCS which has been 
found to be sufficient to reach equilibrium, even in the 
vicinity of the transition. In our scheme, two sets of sim- 
ulations have been performed. In the first one, we pre- 
liminarily explored the thermodynamic behavior of the 
model in order to estimate the position of the maxima 
of the specific heat and susceptibilities and the crossings 
of the fourth order Binder's cumulant. In this first ap- 
proach we used lattice sizes in the interval 20 < L < 50. 
Once the possible transition temperature is determined, 
we refined the results by using single and multiple his- 
togram methods [^l], We produced the histograms 
for each lattice size in the interval 20 < L < 120 and they 
were built at / close to the estimated critical temperatures 
corresponding to the maxima and/or crossing points ob- 
tained in step 1 . Details of the histogram techniques can 
be found in Refs. |l|, ||. 



III. THERMODYNAMIC QUANTITIES AND 
FINITE SIZE SCALING THEORY 



We have devoted our efforts to determine a number of 
thermodynamical quantities, namely, the specific heat, 
magnetization, susceptibility, fourth order Binder's cu- 
mulant and moments of magnetization as described be- 



low. The specific heat is defined as: 

NksT^ ' 



(2) 



where E is the internal energy of the system (computed 
using equation |l|) and N = L?' \s the lattice volume. The 
magnetization is: 



M = — (to) 



where 



\ 




N 



(3) 



(4) 



The susceptibility is defined by the magnetization fluc- 
tuations as: 



(TO 



(TO) 



NkeT 

The fourth order Binder's cumulant reads: 



1 



3(to^ 



(5) 



(6) 



In order to calculate the critical exponent z/, we also de- 
fine the following moments of the magnetization [ p3) : 

(7a) 
(7b) 
(7c) 
(7d) 
(7e) 
(7f) 

where, 
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In critical phenomena the thermodynamic quantities 
are expect to behave in the vicinity of the phase transi- 
tion as pOl pi p5|: 



M ^t^ 



(9a) 
(9b) 
(9c) 
(9d) 



where t = \T — Tc\/Tc is the reduced temperature, M 
is the magnetization, ^ is the correlation length and a, 
/3, 7 and v are critical exponents. Although the criti- 
cal temperature depends on the details of the system in 
consideration, it is observed that the critical exponents 
are universal, depending only on a few fundamental fac- 
tors (20l |23, p5| . The systems are thus divided in a small 
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number of universality classes. Systems belonging to the 
same universality class share the same critical exponents. 
Critical exponents are observed to depend only on the 
spatial dimensionality of the system, the symmetry and 
dimensionality of the order parameter, and the range of 
the interactions within the system. 

In a finite system as those used in Monte Carlo simu- 
lations the divergences in the thermodynamic quantities 
are replaced by smooth functions. Finite size effects are 
therefore of great importance in the analysis of the re- 
sults of Monte Carlo simulations. The theory of finite 
size scaling pO| , p5[ provides one way to extract informa- 
tion concerning the thermodynamic limit properties from 
results obtained in finite systems. The basic assumption 
of this theory is that in the vicinity of the phase tran- 
sition the finite size effects should depend on the ratio 
between the linear dimension of the system and the cor- 
relation length say, L/^ . According to such a theory, 
specific heat, susceptibility and magnetization for a finite 
system, in the vicinity of the phase transition, behave as: 



(10a) 
(10b) 
(10c) 



where M , X and C are proper derivatives of the free en- 
ergy. At Tc {t = 0) these functions are constants and the 
size dependence of specific heat, susceptibility and mag- 
netization follow a pure power law. The size dependence 
of the pseudo-critical temperature, Tc(L), is 



T,iL)^T, + wL-'/'', 



(11) 



where Tc is the critical temperature in the thermody- 
namic limit. Using the size dependence of the magneti- 
zation, equation |l^, and the definition of the moments 
of the magnetization in equation one can easily show 
that such functions behave as: 



Vj w (l/i.)lnL + Vj(<ii/''), 



(12) 



for j = 1,2, ...,6. At t = the functions VjitL^'") are 
constants and then the curves for all Vj have the same 
slope ||2^ providing a very precise method to determine 
both the critical exponent v and the critical temperature. 

Concerning the fourth order Binder's cumulant, it is 
expected that its curves should cross at the same point 
U* = UlT — Tc) for large enough L. Besides that, its 
size dependence is expected to obeys 
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U4 » UiitL^'"). 



IV. RESULTS 



(13) 



In the following we show the results obtained by us- 
ing the histogram method. Each histogram consists of 




FIG. 1: Log-log plot of the maxima of susceptibility as a 
function of the lattice size for L = 20, 40, 80 and 120. The 
error bars are shown inside the symbols. The straight line is 
the best linear fit which gives the exponent 7/1/ = 1.737(1). 




FIG. 2: Specific heat maxima as a function of the lattice size. 
The solid line is the best non-linear fit considering a logarith- 
mic divergence and the dashed line is the best fit considering 
a power law divergence. The error bars are shown inside the 
symbols. 



at least 3 x 10^ configurations. In figure |l| we show a 
log-log plot of the maxima of the susceptibility as a func- 
tion of the lattice size for L = 20,40,80 and 120. The 
data are very well adjusted by a straight line with slope 
7/^ = 1.737(1) exhibiting a power law behavior. The 
specific heat maxima as a function of the lattice size are 
shown in figure ^. In this figure, the solid line represents 
the best non-linear adjust of a logarithmic divergence 
while the dashed one the best power law divergence ad- 
just. It is clear that none of then can adjust our data 
satisfactorily. This result is similar to that obtained for 
the PRM without dipolar interactions, and indicates a 
possible non-divergent specific heat. In figure ^ we show 
the value of 1/v for some temperatures obtained by using 
the moments of the magnetization defined in equation M 
and|l|. Using this method we get tY' = 1.1982(18) and 
1/j/ = 0.74(2). With the value of 1/v we may esti- 
mate the critical temperature using the finite size scaling 
properties of the maxima of the susceptibility and specific 
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FIG. 3: Value of l/i/ obtained by linear fits of Vj versus 
ln(L) for each value of j at different temperatures. Note that 
for T — 1.2000 the value of l/i/ is almost the same for all 
quantities. 
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FIG. 5: Fourth order Binder's cumulant. The critical tem- 
perature was estimated by the crossing point of the largest 
lattice sizes as being Tj^** = 1.203(1). Only a few error bars 
are shown for clarity. 
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FIG. 4: Tc{L) versus L . From a linear adjust we get 
= 1.20022(9) and T,?" = 1.2150(3). 
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FIG. 6: Scaling plot of susceptibility. According to finite 
size scaling theory the susceptibility is expected to 
behave as x ^ L?^'' X(tL ^"'). Note that the curves for dif- 
ferent lattice sizes collapses into a single curve. In the outer 
plot the scaling is done using results from conventional Monte 
Carlo simulations (step 1) for L = 20, 30, 40 and 50. The in- 
set shows the scaling for the histogram results (step 2 in our 
simulations) for L = 20, 40, 80, 120. 



heat, see equation |ri|. In figure ^ we show a plot of Tc{L) 
as a function of L^^". We obtain = 1.20022(9) and 
T^" = 1.2150(3). Using the crossing points of the fourth 
order Binder's cumulant |p6|, see figure ||, we estimate 
the critical temperature T^= 1.203(1). Our best value 
for the critical temperature is thus the mean value of the 
previous estimates discarding the value 

obtained by finite size scaling of the specific heat, since 
its behavior is apparently non-critical. This procedure 
gives Tc = 1.201(1). Plotting In(M^y) versus ln(L) at 
it is possible to obtain the exponent P/v. From a linear 
adjust we get P/i/ = 0.1617(2). In order to verify the 
validity of our results we show in figures ^, (0) and (||), 
the scaling plots of the susceptibility, magnetization and 
fourth order Binder's cumulant according to their finite 
size scaling functions, see equation |l^. Note that all fig- 
ures show a very good collapse of the curves for different 
lattice sizes. 



V. DISCUSSION 

In this work we have studied the phase transition in the 
ferromagnetic dipolar Planar Rotator model (dPRM). 
Our results indicate that the phase transition in this 
model is of order-disorder type and is characterized by 
the exponents ly = 1.277(2), /3 = 0.2065(4) and 7 = 
2.218(5) and by a non-divergent specific heat. Our re- 
sults also indicate that the system have long-range order 
at low temperatures. This conclusion is based on the fol- 
lowing facts: (i) the magnetization for T < Tc does not 
display a significant decrease as the lattice size is aug- 
mented, as for example, has been found in the Rapini's 
work and as expected for a BKT phase transition; 
(m) our results are very well described by a finite size 
scaling theory based on the existence of a low temper- 
ature phase with long-range order and finite correlation 
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FIG. 7; Scaling plot of magnetization. According to finite 
size scaling theory [|o|, ^ this quantity is expected to behave 
as m « L~'^^'"M{tL^^''). Note that the curves for different 
lattice sizes collapses into a single curve. In the outer plot the 
scaling is done using results from conventional Monte Carlo 
simulations (step 1) for L = 20, 30, 40 and 50. The inset 
shows the scaling for the histogram results (step 2 in our 
simulations) for L = 20, 40, 80, 120. 
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FIG. 8: Scaling plot of the fourth order Binder's cumulant. 
According to finite size scaling theory |2^] this quantity is 
expected to behave as U4 ~ UA{tL^^")- Note that the curves 
for different lattice sizes collapses into a single curve. In the 
outer plot the scaling is done using results from conventional 
Monte Carlo simulations (step 1) for L = 20, 30, 40 and 50. 
The inset shows the scaling for the histogram results (step 2 
in our simulations) for L = 20, 40, 80, 120. 



length |2^. In a BKT phase transition there is no 
long-range order in the low temperature phase as a con- 
sequence of the Mermin- Wagner theorem Q. Indeed, the 
results of Maleev |]l5j predict the existence of long-range 
order at low temperatures in the dPRM and our results 
are consistent with this scenario. 

As discussed earlier, recent results by Maier and Schw- 
abl have predicted that this system may belong to 
a new universality class, characterized by the presence 
of long-range order at low temperatures and by an ex- 
ponential behavior of thermodynamic quantities in the 



vicinity of the "critical" temperature. By an exponen- 
tial divergence we mean that the correlation length di- 
verges as the "critical" temperature (Tc) is approached 
as ^ (X exp(6/ {T — Tc)), similarly to the behavior of 
the BKT phase transition, while the behavior of other 
thermodynamic quantities are given by powers of the cor- 
relation length. Nevertheless, our results for the dPRM 
are very well described by power law divergences of ther- 
modynamic quantities. As can be seen in figures |^, |^ 
and we obtained a very good collapse of the curves 
from different lattice sizes for the susceptibility, magne- 
tization and Binder's cumulant. These curves show that 
the critical exponents obtained and the conventional fi- 
nite size scaling theory, that assumes a power law be- 
havior of thermodynamic quantities, describe the Monte 
Carlo data accurately indicating that the phase transi- 
tion in the dPRM is a conventional order-disorder phe- 
nomena with unusual critical exponents. In order to def- 
initely rule out the possibility of this phase transition 
of being in the new universality class proposed by Maier 
and Schwabl, we should make a comparison of our Monte 
Carlo results, using a finite size scaling theory based in 
their predictions, and the conventional finite size scaling 
theory used here. Unfortunately, it is not clear in the 
literature how to obtain a finite size scaling theory for 
exponential divergences. Using a simple replacement of 
the correlation length by the lattice size, which should 
be the first choice, does not give a good collapse of the 
curves, mainly because the determination of the critical 
temperature is quite imprecise in this case and the col- 
lapse of the curves depend appreciably on the value used 
for the critical temperature. In any case, using values for 
the critical temperature close to the maxima of the sus- 
ceptibility we were not able to obtain even a reasonable 
collapse of the curves. 

Once the possibility of this phase transition being in 
the new universality class proposed by Maier and Schw- 
abl is discarded, some questions arise: (i) Why Renor- 
malization Group results do not agree with our Monte 
Carlo simulations? (ii) Is the occurrence of the order- 
disorder transition due to the long-range character of 
dipolar interactions or to some other property of this 
model? A definite answer to these questions may take 
a very long time to be given by virtue of the non-trivial 
characteristics presented by this model. Nevertheless, 
this study gave us some insight about what is happen- 
ing. The RG study of Maier and Schwabl is based 
upon some approximations, for instance, the using of a 
continuous version of dPRM, where the lattice character 
is lost. Since the dipolar interactions have an intrinsic 
anisotropy which depends in a complicated manner on 
the location of each spin in the lattice, the lattice ge- 
ometry could have an strong effect in the system. The 
identification and discussion of the finer points of the RG 
study of the dPRM that cause the discrepancy in the re- 
sults is beyond the scope of this paper. Concerning the 
origin of the order-disorder transition the question is even 
more complicated. The long-range order observed at low 
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temperatures is expected to occur only when full long- 
range interactions are present. Nevertheless, in a recent 
study of the anisotropic Heisenberg model in a bilayer 
system using a cut-off in the dipolar interactions wc 
found the same critical behavior. In fact, the found crit- 
ical exponents {ly = 1.22(9), 7 = 2.1(2) and /3 = 0.18(5)) 
agree inside the errors with those found in this study 
(1/ = 1.277(2), P = 0.2065(4) and 7 = 2.218(5)). This 
observation indicates that the anisotropic character of 
dipolar interactions may be the main responsible by the 
observed critical phenomena. Indeed, this observation 
is not new in the literature. As an example, Fernandez 
and Alonso |^ stated that "Anisotropy has a deeper 
effect on the ordering of systems of classical dipoles in 
2D than the range of dipolar interactions" . In this work 
the authors found that the inclusion of a quadrupolar 
anisotropy changes drastically the phase transition be- 
havior of a system of classical dipoles. Apparently, in 
our system the intrinsic anisotropy of dipolar interac- 
tions play an essential role in the determination of the 
universality class of the dPRM. 

The possible new universality class is not surprising. 
In the theory of critical phenomena pQ, E5| it is ex- 



pected that the critical exponents, and thus the univer- 
sality classes, depend only on the spatial dimensionality 
of the system, the symmetry and dimensionality of the 
order parameter, and the range of the interactions within 
the system, characteristics not shared by the dPRM and 
models of well known universality classes. To the best 
of our knowledge this work and that of Maier and Schw- 
abl are the only ones devoted to the investigation 
of the critical behavior of systems with long-range dipo- 
lar interactions (which are intrinsically anisotropic) and 
exchange interactions in two dimensions. 
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